## Enlistments Data Paper ##
## Prepared by Kevin Fahey ##
## And Doug Atkinson ##
## And Alex Artiles ##
## And Joana Trenaska ##
## Compiled 2023-01-01 Using R 4.2.2

##### PREAMBLE #####

## Clear environment
rm(list=ls())


## Set Seed 
set.seed(20120630)


## Set Working Directory
# setwd() command changed to your working directory


## Load Packages 
library(foreign)
library(tidyverse)
library(broom)
library(vroom)
library(psych)
library(xtable)
library(maptools)
library(usmap)
library(fixest)
library(dotwhisker)


## Read in datasets

# Main dataset for all Figures save Figure 2
dat <- read.csv("2023-01-01 ww2_vols_data.csv")

# Dataset for Figure 2
dat_ww2_women <- read.csv("2023-01-01 dat_ww2_women_all.csv")




##### DESCRIPTIVE STATISTICS #####

## Produce Figure 1 -- Female Volunteer Rates by County, 1943
dat_map <- dat %>%
  dplyr::filter(year == 43) %>%
  mutate(wvol_output = (wvol_prop - mean(wvol_prop, na.rm = T))/sd(wvol_prop, na.rm = T)) %>%
  dplyr::filter(wvol_output <= 3) %>%
  rename(fips = fipscode)
figure1_map <- plot_usmap(regions = "counties", data = dat_map,
                      values = "wvol_output", labels = FALSE,
                      color = "white", size = 0.015,
                      exclude = c("AK", "HI")) +
  scale_fill_continuous(low = "grey20", high = "grey50") +
  labs(title = "",
       fill = "") +
  theme(legend.position = "bottom")
png("figure1_map.png", res = 120)
print(figure1_map)
dev.off()


## Produce Figure 2 -- Female Volunteers, 1941-1945
df_summary_women <- dat_ww2_women %>%
  group_by(enlistyear) %>%
  summarise(a_womenvolunteers = sum(womenvolunteers),
            b_women_whitevolunteers = sum(women_whitevolunteers),
            c_women_blackvolunteers = sum(women_blackvolunteers)) %>%
  pivot_longer(cols = c(a_womenvolunteers, b_women_whitevolunteers, 
                        c_women_blackvolunteers), 
               names_to = "type", 
               values_to = "vols") %>%
  mutate(enlistyear = as.character(paste(enlistyear, sep = ""))) %>%
  filter(enlistyear >= 41 & enlistyear <= 45)
figure2_barchart <- ggplot(df_summary_women,
                                aes(x = enlistyear, y = vols, fill = type)) +
  geom_bar(stat = "identity", position = position_dodge(1.05), width = 0.85) +
  geom_text(aes(label = vols, col = type), 
            vjust = -1.25, size = 2.75, 
            position = position_dodge(0.9)) +
  scale_fill_manual(labels = c("All Women", "White Women", "Black Women"),
                    values = c("grey20", "grey40", "grey60")) +
  scale_color_manual(values = c("grey20", "grey40", "grey60")) +
  guides(color = "none") +
  theme_bw() +
  ylim(0, 60000) +
  labs(x = "Year", y = "Volunteers", fill = "")
png("figure2_barchart.png", res = 100)
print(figure2_barchart)
dev.off()


## Produce Figure 3 -- Female College Graduation Rates by County, 1940
dat_mapcollege <- dat %>%
  rename(fips = fipscode) %>%
  dplyr::filter(year == 43)
figure3_map <- plot_usmap(regions = "counties", data = dat_mapcollege,
                         values = "fgrad_prop", labels = FALSE,
                         color = "white", size = 0.015,
                         exclude = c("AK", "HI")) +
  scale_fill_continuous(low = "grey25", high = "grey55") +
  labs(title = "",
       fill = "") +
  theme(legend.position = "bottom")
png("figure3_map.png", res = 120)
print(figure3_map)
dev.off()


## Produce Figure 4 -- Racial Isolation Rates by County, 1940
dat_mapisolat <- dat %>%
  rename(fips = fipscode) %>%
  dplyr::filter(year == 43)
figure4_map <- plot_usmap(regions = "counties", data = dat_mapisolat,
                        values = "isolation_all", labels = FALSE,
                        color = "white", size = 0.015,
                        exclude = c("AK", "HI")) +
  scale_fill_continuous(low = "grey20", high = "grey50") +
  labs(title = "",
       fill = "") +
  theme(legend.position = "bottom")
png("figure4_map.png", res = 120)
print(figure4_map)
dev.off()


##### REGRESSION OUTPUT #####
model_women <- feols(wvol_prop ~ 
                       mfgwages_prop + retwages_prop + 
                       wholwage_prop + servwage_prop +
                       unemployed_prop + fgrad_prop +
                       totpop40 + urb940_prop +
                       Dem1940 + demcong1940 +
                       dissimilarity_all + isolation_all +
                       wf_prop + wm_prop + negf21_prop + negm21_prop, 
                     data = dat, fixef = c("state", "year"),
                     se = "cluster", cluster = "fipscode")

model_women_white <- feols(wvol_white_prop ~ 
                             mfgwages_prop + retwages_prop + 
                             wholwage_prop + servwage_prop +
                             unemployed_prop + fgrad_prop +
                             totpop40 + urb940_prop +
                             Dem1940 + demcong1940 + 
                             dissimilarity_all + isolation_all +
                             wf_prop + wm_prop + negf21_prop + negm21_prop, 
                           data = dat, fixef = c("state", "year"), 
                           se = "cluster", cluster = "fipscode")

model_women_black <- feols(wvol_black_prop ~ 
                             mfgwages_prop + retwages_prop + 
                             wholwage_prop + servwage_prop +
                             unemployed_prop + fgrad_prop +
                             totpop40 + urb940_prop +
                             Dem1940 + demcong1940 +
                             dissimilarity_all + isolation_all +
                             wf_prop + wm_prop + negf21_prop + negm21_prop, 
                           data = dat, fixef = c("state", "year"), 
                           se = "cluster", cluster = "fipscode")


# Generate tidy output object
models_out <- rbind(
  tidy(model_women) %>% mutate(model = 1),
  tidy(model_women_white) %>% mutate(model = 2),
  tidy(model_women_black) %>% mutate(model = 3)) %>%
  filter(!grepl("factor", term)) %>%
  mutate(model = ifelse(model == 1, "All", 
                        ifelse(model == 2, "White",
                               ifelse(model == 3, "Black", NA))),
         coef_names = rep(c("Manufacturing Wages", "Retail Wages",
                            "Wholesale Wages", "Service Wages", "Unemployment",
                            "Pct. Female Graduates", "Total Population",
                            "Pct. Urban", "Roosevelt Vote Share", 
                            "Dem Cong. Vote Share", "Racial Dissimilarity",
                            "Racial Isolation", "Pct. Women (W)", "Pct. Men (W)",
                            "Pct. Women (B)", "Pct. Men (B)"), 3))


# Separate tidied objects for figures
models_out_fgrad <- models_out %>%
  filter(term == "fgrad_prop") %>%
  mutate(term = rep(c("Pct Female Graduates"), 3))
models_out_race <- models_out %>%
  filter(term == "dissimilarity_all" | term == "isolation_all") %>%
  mutate(term = rep(c("Racial Dissimilarity", "Racial Isolation"), 3))
models_out_econ <- models_out %>%
  filter(term == "mfgwages_prop" | term == "retwages_prop" | term == "wholwage_prop" | 
           term == "servwage_prop" | term == "unemployed_prop") %>%
  mutate(term = rep(c("Manufacturing Wages", "Retail Wages", "Wholesale Wages",
                      "Service Wages", "Unemployment"), 3))
models_out_demo <- models_out %>%
  filter(term == "wf_prop" | term == "wm_prop" | 
           term == "negf21_prop" | term == "negm21_prop" |
           term == "urb940_prop") %>%
  mutate(term = rep(c("Pct. Urban", "Pct. Women (W)", "Pct. Men (W)",
                      "Pct. Women (B)", "Pct. Men (B)"), 3))
models_out_poli <- models_out %>%
  filter(term == "Dem1940" | term == "demcong1940") %>%
  mutate(term = rep(c("Roosevelt Vote Share",
                      "Dem Cong. Vote Share"), 3))


## Produce Figure 5 -- Coefficient Plot of Influence of College Education on Female Volunteerism, with 95% CIs
figure5_dwplot <- dwplot(x = models_out_fgrad, dodge_size = 0.5, ci = 0.95,
                              dot_args = list(size = 3),
                              whisker_args = list(size = 1, aes(linetype = model)),
                              vline = geom_vline(xintercept = 0.0,
                                                 colour = "black",
                                                 linetype = 1)) +
  scale_colour_manual(labels = c("Black Women", "White Women", "All Women"),
                      values = c("grey60", "grey40", "grey20")) +
  theme_bw() + xlab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.y = element_text(angle = 45, size = 12),
        axis.text.x = element_text(angle = 25, size = 12),
        plot.margin = margin(1, 1, 1, 1, "cm")) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom")
png("figure5_dwplot.png", res = 95)
print(figure5_dwplot)
dev.off()


## Produce Figure 6 -- Coefficient Plot of Influence of Racial Variables on Female Volunteerism, with 95% CIs
figure6_dwplot <- dwplot(x = models_out_race, dodge_size = 0.5, ci = 0.95,
                             dot_args = list(size = 3),
                             whisker_args = list(size = 1, aes(linetype = model)),
                             vline = geom_vline(xintercept = 0.0,
                                                colour = "black",
                                                linetype = 2)) +
  scale_colour_manual(labels = c("Black Women", "White Women", "All Women"),
                      values = c("grey60", "grey40", "grey20")) +
  theme_bw() + xlab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.y = element_text(angle = 45, size = 12),
        axis.text.x = element_text(angle = 25, size = 12),
        plot.margin = margin(1, 1, 1, 1, "cm")) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom")
png("figure6_dwplot.png", res = 95)
print(figure6_dwplot)
dev.off()


## Produce Figure 7 -- Coefficient Plot of Influence of Economic Variables on Female Volunteerism, with 95% CIs
figure7_dwplot <- dwplot(x = models_out_econ, dodge_size = 0.5, ci = 0.95,
                             dot_args = list(size = 3),
                             whisker_args = list(size = 1, aes(linetype = model)),
                             vline = geom_vline(xintercept = 0.0,
                                                colour = "black",
                                                linetype = 2)) +
  scale_colour_manual(labels = c("Black Women", "White Women", "All Women"),
                      values = c("grey60", "grey40", "grey20")) +
  theme_bw() + xlab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.y = element_text(angle = 45, size = 12),
        axis.text.x = element_text(angle = 25, size = 12),
        plot.margin = margin(1, 1, 1, 1, "cm")) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom")
png("figure7_dwplot.png", res = 95)
print(figure7_dwplot)
dev.off()


## Produce Figure 8 -- Coefficient Plot of Influence of Political Variables on Female Volunteerism, with 95% CIs
figure8_dwplot <- dwplot(x = models_out_poli, dodge_size = 0.5, ci = 0.95,
                             dot_args = list(size = 3),
                             whisker_args = list(size = 1, aes(linetype = model)),
                             vline = geom_vline(xintercept = 0.0,
                                                colour = "black",
                                                linetype = 2)) +
  scale_colour_manual(labels = c("Black Women", "White Women", "All Women"),
                      values = c("grey60", "grey40", "grey20")) +
  theme_bw() + xlab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.y = element_text(angle = 45, size = 12),
        axis.text.x = element_text(angle = 25, size = 12),
        plot.margin = margin(1, 1, 1, 1, "cm")) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom")
png("figure8_dwplot.png", res = 95)
print(figure8_dwplot)
dev.off()


## Produce Figure 9 -- Coefficient Plot of Influence of Demographic Variables on Female Volunteerism, with 95% CIs
figure9_dwplot <- dwplot(x = models_out_demo, dodge_size = 0.5, ci = 0.95,
                             dot_args = list(size = 3),
                             whisker_args = list(size = 1, aes(linetype = model)),
                             vline = geom_vline(xintercept = 0.0,
                                                color = "black",
                                                linetype = 2)) +
  scale_colour_manual(labels = c("Black Women", "White Women", "All Women"),
                      values = c("grey60", "grey40", "grey20")) +
  theme_bw() + xlab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.y = element_text(angle = 45, size = 12),
        axis.text.x = element_text(angle = 25, size = 12),
        plot.margin = margin(1, 1, 1, 1, "cm")) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom")
png("figure9_dwplot.png", res = 95)
print(figure9_dwplot)
dev.off()
